#Joseph Simons
#started on 6.29.2014

#######################################################
#This script was most recently modified on 3.15.2015
#######################################################

##This script runs replications of the shotts models
#using the 2012 Simons-Mallinson data and the corrected 
#script for aggregating the data and merging it with the
#DW-Noms data. 

#This script also runs the original and modified interactive
#specifications of the original simons-malinson model.


#Last modified on 2.15.2015 at 2024 hrs.


#################


rm(list=ls())
#wd<-  #set the working directory where you intend to store the replication data
#setwd(wd)  


jopnoms.south<-read.csv("simons.malinsonreplicationdata.csv") #call up the data


#t(10%)=1.65
#t(5%)=1.95
#note that R computes the one tailed p-value for a given alpha
#if you want the two tailed value, you need to divide the alpha by 2
#eg, if you want a two-tailed test with a 0.05 rejection region, you need use 
#a value equal to 0.05/2=0.025
#for example, p=0.95 is really giving you the 90% critical value for the 
#tailed test.

#Marginal Effects and variance
mfx.fmmdists<-function(I.MS,repgerr,bipartgerr,alpha,hccm=FALSE){
	coeffs<-coefficients(I.MS)
	if(hccm==TRUE){require(car)
		vcv<-hccm(I.MS)
	}else{
		vcv<-vcov(I.MS)
	}
	mfx.fmmdists<-coeffs[2]+coeffs[9]*(repgerr)+coeffs[10]*(bipartgerr)
	var.mfx.fmmdists<-vcv[2,2]+(repgerr^2)*vcv[9,9]+(bipartgerr^2)*vcv[10,10]
		+2*(repgerr)*vcv[2,9]+2*(bipartgerr)*vcv[2,10]+2*(repgerr*bipartgerr)*vcv[9,10]
	se.mfx<-sqrt(var.mfx.fmmdists)
	t<-mfx.fmmdists/se.mfx
	crit<-qt(df=I.MS$result$df.residual,p=alpha)
	stats<-c(mfx.fmmdists,var.mfx.fmmdists,se.mfx,t,crit)
return(stats)
}


f.MS2<-as.formula(flibs.Shotts~fmmdists+repgerr+bipartgerr+votlib+blackpop+hispop
	+flibs.Shotts.lag+I(fmmdists*repgerr)+I(fmmdists*bipartgerr))

f.shotts<-as.formula(flibs.Shotts~fmmdists+demgerr+votlib
	+blackpop+hispop+flibs.Shotts.lag)

#####################################
####Run the models###################
#####################################

Z.MS2.110<-zelig(formula=f.MS2,model="normal",
	data=jopnoms.south, weights=jopnoms.south$districts)

Z.shotts.110<-zelig(formula=f.shotts,model="normal",
	data=jopnoms.south, weights=jopnoms.south$districts)

summary(Z.MS2.110)
mfx.fmmdists(Z.MS2.110,repgerr=1,bipartgerr=0,alpha=.05) #repgerr
mfx.fmmdists(Z.MS2.110,repgerr=0,bipartgerr=0,alpha=.05) #demgerr
mfx.fmmdists(Z.MS2.110,repgerr=0,bipartgerr=1,alpha=.05)#bipartgerr

summary(Z.shotts.110)


##Replicate the oringinal results from
#shotts 2003 using the dataset shotts.data.south

Z.shotts.50<-zelig(formula=f.shotts,model="normal",
	data=shotts.data.south, weights=shotts.data.south$districts)

#Estimate the original shotts model using the original data#####
################################################################
Z.shotts.50.2003<-zelig(formula=f.shotts,model="normal",
	data=withlags.south, weights=withlags.south$Districts)

#############################################
#make the graphics for the interactive model#
#############################################


fmmdists.range<-seq(min(fmmdists),max(fmmdists),by=0.01)
fmmdists.shift<-mean(fmmdists)+sd(fmmdists)


par(mfrow=c(1,3))
xc.range1<- setx(Z.MS2.110,repgerr=1,bipartgerr=0, 
	fmmdists=fmmdists.range,
	flibs.Shotts.lag=mean(jopnoms.south$flibs.Shotts.lag,na.rm=T))    #a republican gerrymander
xc.range2<- setx(Z.MS2.110,repgerr=0,bipartgerr=1, 
	fmmdists=fmmdists.range,
	flibs.Shotts.lag=mean(jopnoms.south$flibs.Shotts.lag,na.rm=T))    #a bipartisan gerrymander
xc.range3<- setx(Z.MS2.110,repgerr=0,bipartgerr=0, 
	fmmdists=fmmdists.range,
	flibs.Shotts.lag=mean(jopnoms.south$flibs.Shotts.lag,na.rm=T))    #a democrat gerrymander gerrymander
s.out1 <-sim(Z.MS2.110, x=xc.range1)
s.out2 <-sim(Z.MS2.110, x=xc.range2)
s.out3 <-sim(Z.MS2.110, x=xc.range3)

ylab<-"Expected Values (flibs|X)"
xlab<-"Range of fmmdists"

setEPS()
postscript("repev110.eps")
plot.ci(s.out1,qi="ev",main="Republican Gerrymander",ylab=ylab, xlab=xlab,var="fmmdists",xlim=c(0,0.3))
#dev.copy2pdf(file="repev110.pdf")
dev.off()

setEPS()
postscript("bipartev110.eps")
plot.ci(s.out2,qi="ev",main="Bipartisan Gerrymander",ylab=ylab, xlab=xlab,var="fmmdists",xlim=c(0,0.3))
#dev.copy2pdf(file="bipartev110.pdf")
dev.off()

setEPS()
postscript("demev110.eps")
plot.ci(s.out3,qi="ev",main="Democrat Gerrymander",ylab=ylab, xlab=xlab,var="fmmdists",xlim=c(0,0.3))
#dev.copy2pdf(file="demev110.pdf")
dev.off()
